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ABSTRACT 


The  study  of  lateral  variations  of  earth  structure  has  been  stimulated  recently 
by  several  factors,  especially  the  theory  of  plate  tectonics  and  the  increasing  use  of 
large  seismic  arrays.  The  extension  of  seismic  ray  theory  to  two  and  three  dimen¬ 
sional  structures  is  thus  of  great  practical  importance. 

The  problem  of  ray  tracing  in  a  generally  heterogeneous  medium  is  treated,  using 
the  calculus  of  variations  and  Fermat's  principle  of  stationary  time.  The  solution  is 
expressed  in  terms  of  a  system  of  five  simultaneous  first  order  differential  equations 
giving  the  variation  with  time  of  the  position  and  direction  of  motion  of  a  point  on  a 
ray  in  terms  of  the  wave  speed  and  its  spatial  derivatives  in  the  medium.  The  form  of 
the  equations  is  particularly  convenient  for  solution  on  a  digital  computer. 

The  effect  upon  the  wave  amplitudes  of  geometric  spreading  can  also  be  calculated 
by  the  inclusion  of  ten  additional  equations  in  the  system,  although  this  greatly  increases 
the  computational  labor. 

If  the  earth  model  has  certain  symmetry  properties,  then  constants  of  the  motion 
along  each  ray  can  be  found  which  simplify  the  calculations.  For  example,  for  two- 
dimensional  models  in  which  the  velocity  depends  only  on  the  coordinates  r  (radius) 
and  0  (colatitude)  one  equation  can  be  eliminated  from  the  ray  tracing  system,  and 
two  more  can  be  eliminated  from  the  amplitude  equations. 

The  propagation  of  surface  waves  on  an  earth  with  geographical  variations  can 
be  treated  by  a  simplified  special  case  of  the  method  presented  here. 

Accepted  for  the  Air  Force 

Joseph  R.  Waterman,  Lt.  Col.,  USAF 

Chief,  Lincoln  Laboratory  Project  Office 
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Introduction 


In  the  past,  most  studies  of  seismic  waves  have  been  based  upon  the  assump¬ 
tion  that  the  earth  is  spherically  symmetrical,  with  physical  properties  depending 
only  on  radius.  Recently,  however,  it  has  become  apparent  that  this  approximation 
is  often  not  justified  Oulian,  1969)  and  more  attention  has  been  directed  by  seismologists 
to  the  study  of  lateral  variations  of  structure.  Two  developments  in  particular  have 
stimulated  interest  in  wave  propagation  through  two  and  three  dimensional  structures: 
the  existence  of  large  seismic  arrays,  which  have  yielded  clear  evidence  of  azimuthally 
dependent  travel  time  anomalies'*',  and  the  theory  of  plate  tectonics  which  hypothesizes 
the  existence  of  relatively  systematic  large  scale  motions  in  the  upper  part  of  the 
earth  (Isacks  et.  al.  ,  1968;  Davies  and  McKenzie,  1969).  This  paper  presents  a  for¬ 
mulation  of  seismic  ray  theory  which  is  applicable  to  such  problems.  The  system  of 
differential  equations  derived  is  easily  soluable  numerically  to  trace  rays  through 
complicated  structures. 

Ray  Tracing 

Recently  V.  A.  Eliseevnin  (1965)  has  formulated  the  ray  problem  for  an  arbi¬ 
trarily  inhomogeneous  medium.  Starting  with  the  Eikonal  equation,  he  derived  the 
following  system  of  six  simultaneous  differential  equations  for  the  motion  of  a  distur¬ 
bance  along  a  ray: 

*Niazi,  1966;  Bolt  and  Nuttli,  1966;  Nuttli  and  Bolt,  1969;  Greenfield  and  Sheppard,  1969. 
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where: 


x,y,z  are  the  cartesian  coordinates  of  the  distrrbance  at  a  particular 
time. 


a,  g,  y  are  the  direction  angles  of  the  tangent  to  the  ray. 


v(x,y,z)  is  the  wave  speed. 

dv  dv  dv 

The  partial  derivatives  ^  ^  ,  and  ^  indicate  spatial  derivatives 

atx,  y,  z,  not  derivatives  along  the  ray  path. 

Only  five  of  these  equations  are  independent,  because  the  angles  a,  8  and  y  are 


connected  by  the  relation  cos2  a  +  cos2  8  +  cos2  y  =  1.  A  much  simpler  way  of  stating 
these  equations  follows  if  we  define  a  slowness  vector,  S,  such  that 
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=  v2S 
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The  components  of  S  are 
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The  rate  of  change  of  S  is 


-  (S-V  v)  v  S 


where  e  ,  e  and  e 
x  y  z 

(1)  this  becomes 


are  unit  vectors  in  the  direction  of  the  coordinate  axes  and,  using 


S 


(lb) 


Despite  the  simplicity  of  (la)  and  (lb),  for  computational  purposes  expressions 
such  as  (1)  in  terms  of  angles  instead  of  slowness  are  usually  more  useful,  particularly 
in  curvilinear  coordinate  systems.  We  shall  give  a  different  derivation  of  equations 
(1),  based  on  Fermat's  principle  of  least  time,  and  carry  out  the  derivation  in  spheri¬ 
cal  coordinates,  so  that  the  result  will  be  in  a  seismologically  useful  form. 

Let  r,  0,  0,  be  the  spherical  coordinates,  at  time  t,  of  a  point  on  a  ray.  Further, 

letting  §  ,  e.,@  ,  be  the  conventional  unit  vectors  for  spherical  coordinates,  define: 

r  9  0 

i  =  angle  between  ray  direction  and  e  . 

i  =  angle  between  ray  direction  and  ea. 

9  9 

i  =  angle  between  ray  direction  and  e  . . 

0  0 

The  first  three  differential  equations  follow  geometrically: 
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To  find  the  differential  equations  for  i  ,  ifl ,  and  i  ,  we  consider  conditions  for  the 
travel  time  to  be  stationary  with  respect  to  small  changes  in  the  ray  path.  Using  9 
as  a  parameter,  let  the  ray  path,  c,  be  specified  in  terms  of  r(0)  and  0(9).  The 
travel  time  of  a  ray  between  two  points  where  9  =  0!  and  9  =  0S  is 

T  =  r  la  _  p  r  d9 

J  v  J  v  cos  i„ 
c  A  0 


where  do  is  an  element  of  length  of  the  ray  path.  Consider  a  small  change  in  the  ray 
path  specified  by  6r(0),  6  0(9)  with  6r(9].)  =  6r(92)  =  6  0 (0 x)  =  6  0(02)  =  0,  that  is,  with 
the  end  points  of  the  ray  fixed.  The  change  in  the  travel  time  is 
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Performing  some  algebraic  manipulations,  and  integrating  by  parts,  this  can  be 
written 
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As  before,  and  are  ordinary  derivatives  with  respect  to  position,  not 

derivatives  along  the  ray.  Since  we  want  6T  =  0  for  arbitrary  6r  and  60,  the  coef¬ 
ficients  of  6r  and  60  in  (6)  must  vanish.  Noting,  from  (2),  (3),  and  (4)  that 


dr  _  r  cos  ir 
d0  cos  ig 

and 


.  ,  cos  1 
d0  _  0 
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these  conditions  can  be  written  as 
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Using  (3),  we  see  that  these  are  the  expressions  for  i  and  i  ,  .  The  expression  for 

r  0 

i„  can  be  found  from  them  using  the  fact  that  cos3i  +  cos3i.  +  cos3i  =  1. 

0  6  r  0  0 

Instead  of  i  ,  or  iQ ,  it  is  simpler  to  use  the  angle  G  between  the  vertical  plane 

0  C7 
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tangent  to  the  ray  and  the  meridional  plane.  We  have 


cos  1 

cos  C  =  — — — 
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cos  1 
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sin  i 


so 


.  sin  in  d  in  cos  i0cos  i  d  i 

9  0  ^  0  r  r 
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r  sin  i  9  0  sin  i 
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- : - ~  ~  am  i  t 
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and  we  can  write  the  five  equations  for  the  ray  as 
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Amplitudes  -  Geometric  Spreading 

Two  phenomena  affect  the  amplitudes  of  body  waves:  geometric  spreading  of 
the  rays  and  attenuation  due  to  anelasticity.  We  will  direct  our  attention  to  geometric 
spreading  first,  assuming  the  earth  is  non-dissipative. 
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Let 


1(1^  Cq)  =  power/nnit  solid  angle  radiated  at  the  focus 
E(®,  $)  =  power/unit  area  of  wavefront  at  the  point  of  observation 

where  if0  CQ  are  the  initial  values  of  if,  and  ©,  §  are  the  values  of  0 , 0  at  the  point 
of  observation.  In  a  non -dissipative  earth 


I(iro>C0)dfi  =  E©,§)dA  (14) 

where  dQ  and  dA  are  the  corresponding  elements  of  solid  angle  at  the  source  and  sur¬ 
face  area  of  wavefront  at  the  receiver  and  are  given  by 


dQ  =  sin  i  di  dC 
ro  ro  o 

. .  _  R2sin© 

dA  = - r=-  d©d§ 

cos  1 


(15) 


R  is  the  earth's  radius.  Here  i  refers  to  the  value  at  the  observation  point.  d©and 
d$  refer  to  changes  with  r  held  fixed,  along  the  earth's  surface,  not  along  a  wavefront, 
d  i  ,  dCQ,  d©  and  d$  are  related  by  the  Jacobian  of  the  transformation  from  ©,  f  to 
*ro'  ^o  ^e^ed  by  die  rays: 
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From  (14),  (15),  and  (16)  we  get 
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To  evaluate  the  partial  derivatives  in  (16)  we  must  solve  ten  more  differential  equations, 
dr  dr 


for 


di 


,  d  »  •  •  •  ,  -=tt£  simultaneously  with  (91-(13).  These  equations  are 

bo  dlro  *o 

obtained  by  differentiating  equations  (9)-(13)  with  respect  to  iro  and  CQ  and  reversing 

the  order  of  differentiation  -  [4£]  =  iL  ],  etc. ),  yielding 
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and  we  have  used  q  to  represent  either  of  the  initial  angles  irQ  or  The  derivatives 

tt —  ,  .  -■■■  .  etc.  thus  obtained  are  those  which  apply  when  the  travel  time  is  held 
ro  ^o 

fixed;  that  is  they  apply  to  values  on  a  particular  wavefront.  Derivatives  along  the 
earth's  surface  (r  constant)  needed  in  equation  (16)  can  be  obtained  using  the  identity 


JLl  =  JL  I  -  (If  I  /  l!  I  )  I 

Sq  'r  Sq 't  ^Sq  't  '  St  'q;  St  1  q  . 


(24) 


The  constant  q  derivatives  are  just  those  given  in  equations  (9)-(13). 
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Amplitudes  -  Attenuation 


The  attenuation  of  seismic  waves  due  to  anelasticity  is  quite  easy  to 
calculate  if  the  intrinsic  quality  factor,  Q,  is  known  as  a  function  of  position  and  the 
ray  path  has  been  calculated.  The  power  in  the  wave  is  simply  reduced  by  the  factor 


where  oi  is  the  angular  frequency  of  the  wave  and  the  integral  is  evaluated  along  the 
ray  path.  The  power  per  unit  area  of  wavefront,  P,  is  related  to  the  amplitude,  A, 


by 


P  =  ■§  p  Voj  s  A2 


for  both  compressional  and  shear  waves,  where  p  is  the  density  of  the  medium. 


Symmetry  relations 

In  special  cases  in  which  the  earth  model  posesses  certain  symmetry  properties, 
the  5  ray  variables  are  no  longer  independent  and  one  or  more  of  the  equations  (9) -(13) 
may  be  eliminated  from  the  system  to  be  solved.  For  example,  in  cartesian  coordinates, 
if  the  velocity  is  independent  of  one  of  the  coordinates  then  by  (lb)  the  corresponding 
component  of  the  slowness  vector  S  is  constant  along  a  ray.  In  spherical  coordinates 
the  situation  is  core  complicated.  The  rate  of  change  of  S  is 


S  =  [Sr-SQe-S0sin  00]  er 


+  [Sq+S^.0  -S^cos  0  0]  eft 
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and,  with  (lb)  this  gives 
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From  (la),  (9),  (10),  and  (11)  the  components  of  S  are 


cos  i 
_ r 

v 

sin  i  cos  C 
r 

v 

sin  i  sin  C 
r  * 

v 


These  equations  can  be  used  to  verify  the  following  relations : 

V  =  f(r):  (rS0)2  +  (rS^)3  =  (£  sin  ir)3  =  const.  (25) 

v  =  f(r,  0):  r  sin  0  S  ^,  =  ~  sin  i^  sin  Q  sin  0  =  const.  (26) 

The  amplitude  calculations,  also,  can  be  simplified  by  using  symmetry  rela¬ 
tions.  For  example,  when  the  velocity  depends  only  on  r  and  0,  from  equation  (26) 
the  quantity 

p  =  —  sin  i  sin  C  sin  0 
v  v  r  * 
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is  constant  along  a  ray,  and  thus  so  are  its  derivatives 


dp 

dq 


_  rl  dr  1  Dv  ,  „  .  hlr  ,  „  -  dC  ^  „  A  d  0  , 

^  Lr  dq  v  Dq  r  dq  *  dq  dq 


(27) 


Dv 

where,  as  before,  q  represents  either  irQ  or  and  ^  is  defined  in  (23).  It  further 
follows  that  the  quantity  in  brackets  in  (27)  is  constant  and  this  fact  can  be  used  to 
eliminate  two  of  the  equations  (18) -(22)  from  the  amplitude  calculations. 


Ray  Tracing  for  Surface  Waves 

It  should  be  pointed  out  that  a  very  similar  approach  can  be  used  to  calculate 
surface  wave  paths  on  an  earth  model  in  which  the  surface  wave  velocity  is  a  function 
of  geographic  position.  The  ray  tracing  equations  are  obtained  from  equations  (10), 
(11)  and  (13)  by  setting  r  =  R  and  i  =  TT/2: 


cos  £ 

(28) 

v  .  - 

R  sin  0 

(29) 

sin  C  dv  _  cos  Q 

R  d0  R  sin  0 

dv  v  .  -  _  _ 

50  •  R  SU1  i  COt  9 

(30) 

The  calculation  of  amplitudes,  too,  follows  a  development  similar  to  that 
given  above.  By  analogy  with  equation  (14)  we  have 

I(CQ)dC0  =  E  (0,0)  dl  (31) 

where  I  (£  )  =  power/unit  angle  radiated  at  the  source 

E(0,0)  =  power/unit  length  of  wavefront,  dl,  at  0,0. 
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Using  the  fact  that 


hi 

equation  (31)  gives 

E  = 


rJ  (|£  )2  +  sin3  e  (||  y 


and  the  additional  differential  equations  needed  to  calculate  and  ^  can  be  ob¬ 
tained  from  equations  (28)-(30)  the  same  way  equations  (18)-(22)  were  derived  from 
<9)-<13). 


Examples 

As  examples  of  the  use  of  the  ray  tracing  technique,  we  will  show  ray  paths 
calculated  from  two-dimensional  models  of  dipping  high-velocity  slabs  within  the 
earth's  upper  mantle,  such  as  are  present  beneath  island  arcs.  The  numerical  solu¬ 
tion  of  the  differential  equations  has  been  carried  out  using  the  step  size  extrapola¬ 
tion  method  of  Bulirsch  and  Stoer  (1966),  which  is  particularly  fast  and  accurate. 

Figure  1  shows  ray  paths  and  wavefronts  for  a  very  simple  model,  in  which 
the  velocity  is  specified  by  an  analytical  function 

v  =  vq  +  A  exp  [-£)s  -  j  ] 

where  x  is  the  distance  from  the  slab  axis  (assumed  to  dip  at  a  45°  angle,  and  z  is 
depth  from  the  earth’s  surface,  v  ,  A,  h,  and  d  are  constants.  In  this  example,  the 
earthquake  focus  is  located  near  the  edge  of  the  slab,  thus  producing,  in  addition  to  a 
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Fig.  1 .  P  wave  ray  paths  and  wavefronts  calculated  for  a  focus  near  the  edge 
of  a  simple  analytical  model  of  a  hi gh -velocity  lithospheric  slab  (see  text). 
The  slab  parameters  are:  dip  =  45°;  velocity  outside  slab,  vQ  =  8.  0  km/sec; 
maximum  velocity  contrast,  A  =  0.  8  km/ sec;  slab  thickness,  2  h  =  80  km; 
depth  of  penetration,  d  =  300  km. 
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large  shadow  zone,  a  region  of  strong  focusing  of  energy. 

Figure  2  shows  ray  paths  and  wavefronts  for  a  wave  emerging  from  the  mantle 
beneath  a  more  detailed  slab  model.  The  velocity  model  is  derived  from  the  theore¬ 
tical  temperature  field  calculated  numerically  by  Toksoz  et  al  (1971).  In  this  case, 
the  velocity  is  specified  at  a  grid  of  points,  and  interpolation  between  these  points  is 
accomplished  using  cubic  splines,  so  that  the  velocity  and  its  first  spatial  derivatives 
are  continuous.  Again,  a  shallow  zone  is  produced,  along  with  regions  of  focusing. 
Seismic  stations  on  or  near  island  arcs  would  thus  be  expected  to  have  "blind  spots" 
for  earthquakes  in  certain  regions. 

Conclusions 

The  formulation  of  ray  theory  given  here  is  applicable  to  a  medium  with  any 
type  of  inhomogeniety,  and  is  well  suited  for  numerical  computations.  It  allows  the 
computation  of  amplitudes,  considering  both  geometric  and  anelastic  effects,  though 
the  computational  labor  may  be  greatly  increased.  Furthermore,  it  can  be  extended 
easily  to  problems  such  as  surface  wave  propagation  on  an  earth  with  geographical 
variations. 
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Fig.  2.  Ray  paths  and  wavefronts  for  P  wave  emerging  beneath  an  island  arc. 
The  slab  model  is  derived  from  numerical  temperature  field  calculations. 
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